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^ ■ We propose an efficient, accurate method to integrate the basins of attraction of 

a smooth function defined on a general discrete grid, and apply it to the Bader 
charge partitioning for the electron charge density. Starting with the evolution of 
trajectories in space following the gradient of charge density, we derive an expression 
for the fraction of space neighboring each grid point that flows to its neighbors. 
This serves as the basis to compute the fraction of each grid volume that belongs 
to a basin (Bader volume), and as a weight for the discrete integration of functions 
over the Bader volume. Compared with other grid-based algorithms, our approach is 
robust, more computationally efficient with linear computational effort, accurate, and 
has quadratic convergence. Moreover, it is straightforward to extend to non-uniform 
>• ■ grids, such as from a mesh-refinement approach, and can be used to both identify 

basins of attraction of fixed points and integrate functions over the basins. 
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I. INTRODUCTION 



Based on density functional theory (DFT)A calculations, decomposing the charge or the 
energy of a material into contributions from individual atoms can provide new information 
for material properties. Bader's "atoms in molecules" theory provides an example of a 
partitioning based on the charge density, and following the gradient at a particular point 
in space to the location of a charge density maximum centered at an atom — defining basins 
of attraction of fixed points of the charge density. Bader defines the atomic charges and 
well-defined kinetic energies as integrals over these Bader volumes, 2 Q p . Each Bader volume 
contains a single electron density maximum, and is separated from other volumes by a 
zero flux surface of the gradients of the electron density, Vp(r) • h = 0. Here, p(f) is the 
electron density, and n is the unit vector perpendicular to the dividing surface at any surface 
point f e dVtp. Each volume Q p is defined by a set of points where following a trajectory 
of maximizing p reaches the same unique maximum (fixed point). In practical numerical 
calculations, where the charge density is defined on a discrete grid of points in real space, it 
is very challenging to have an accurate determination of a zero flux surface. 

Different approaches for condensed, periodic systems have relied on analytic expressions 
of the density^ or discretizing the charge density trajectories^-. Early algorithms were 
based on the electron density calculated from analytical wavefunctions of small molecules, 
and integration along the gradient paths. Most current developments are based on a grid of 
electron density, which is important for DFT calculation and also applicable to analytical 
density function of small molecules. One octal tree algorithm^ uses a recursive cube sub- 
division to find the atomic basins robustly, but practically is not applicable to complicated 
topologies due to huge computational cost. The "elastic sheet" method 6 defines a series of 
fictitious particles which gives a discrete representation of zero-flux surface. Particles are re- 
laxed according to the gradients of charge density and interparticle forces. This method will 
not work for complex surface with sharp cusps or points. Recently, Henkelman et al. devel- 
oped an on-grid method^ to divide an electron density grid into Bader volumes. This method 
can be applied to the DFT calculations of large molecules or materials. They discretize the 
trajectory to lie on the grid, ending at the local maximal point of the electron density. The 
points along each trajectory are assigned to the atom closest to the end point. Although this 
method is robust, and scales linearly with the grid size, it introduces a lattice bias caused by 
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the fact that ascent trajectories are constrained to the grid points. The near-grid method 8 
improves this by accumulating a correction vector — the difference between the discretized 
trajectory and the true trajectory — at each step. When the correction vector is sufficiently 
large, the discrete trajectory is corrected to a neighboring grid point. This method corrects 
the lattice bias, and also scales linearly with respect to the size of grids. However, both 
grid trajectory methods require iteration to self-consistency in volume assignments. Also, 
the integration error scales linearly with the grid spacing, so very fine grids are required 
in numerical calculations to provide the correct Bader volume, reducing its applicability for 
accurate calculations in a large system. Lastly, a new algorithm uses a "divide and conquer" 
adaptive approach with tetrahedra; tetrahedra are continuously divided at the boundaries 
of Bader volumes, with the weight of each tetrahedra given by the number of vertices that 
belong to each volume 9 . Such an approach retains linear scaling with the grid spacing, but 
requires mesh refinement near boundaries to deal with the linear convergence of the error 
with the grid spacing. 
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FIG. 1. A schematic illustration of the zero flux surface, the near-grid algorithm, and weighted 
integration. The zero flux dividing surface separate volumes A and B, where arrows denote charge 
density gradients (a). The normal component is zero for any point on the surface V/3 • ft = 0. The 
near-grid method 8 gives grid-based partition (b), however energy density integration based on this 
grid-based partition would cause integration error due to finite grid sizes. A weight function (c) 
representing volume fractions of the cell of each grid point is introduced to reduce the error due to 
a finite grid. 

Figure [JJ illustrates the partition of the real space into two Bader volumes A and B by 
a zero flux dividing surface. The component of Vp along the surface normal n is zero for 
any point on the surface dA or OB. A grid-based partition algorithm, such as the near-grid 
method, divides space into volume surrounding around each grid point, and assigns each 
grid volume to a particular Bader volume. Even though grid points may be assigned to 
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Bader volumes correctly, the density integration based on the grid-based partition would 
bring in numerical integration error that scales linearly with the grid spacing. Introducing 
a "weight" integrand representing the fraction of grid volume that belongs to a particular 
Bader volume smooths out the grid-based partition, and improves the integration accuracy 
and scales quadratically with the grid spacing. The atomic contribution is neither 1 nor 
at the dividing surfaces, but fractional. In Figure HJ red represents a weight of 1 to atom A 
for grid points closer to atom A, and transitions to white for a weight of for grid points 
away from atom A. 

The grid-based weight representing volume fractions of each grid volume assigned to 
different atoms, and gives a more accurate integrand for the integration of either charge 
density or of kinetic energy over the Bader volumes. The weight is computed from the 
total integrated flux of trajectories in a grid volume to neighboring grid volumes. The al- 
gorithm is robust, efficient with linear computing time in the number of grid points, and 
more accurate than other grid-based algorithm. Surprisingly, it combines both better error 
scaling — quadratic in the grid spacing — and improved computational efficiency. Moreover, 
it is straightforward to apply to nonuniform grids, such as would result from an adap- 
tive mesh-refinement approach. In Section [Til we derive the algorithm to constructing a 
grid-based weight to perform numerical integrals over Bader volumes. Section ITTTI presents 
examples including three dimensional charge density from three Gaussian functions in FCC 
cell, Ti02 bulk, and NaCl crystal. Finally, we show the improved computational efficiency 
in Section IIVI The end result is a simple, extendable, computationally efficient algorithm 
with quadratic integration error. 

II. THE WEIGHT METHOD 

The Bader partitioning of space defines volumes by the endpoint of a trajectory following 
the gradient flow of the charge density, Vp. We assume that p has continuous first and 
second derivatives throughout all space of interest, and has a set of discrete local maxima 
(fixed points) xx, %2, etc., where Vp = and the matrix VVp is negative-definite. The basin 
of attraction A n , of a fixed point x n is the set of points which flow to the fixed point x n along 
the charge density gradient. That is, for any point r, we can integrate the trajectory given 
by x(t) = Vp(af), with the initial condition x(0) = r, to find lim^oo x(t). Each trajectory 
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will end at fixed point x n , and except for a set of points with zero volume in space, the 
extremum is a local maximum; the basin of attraction A n are all points r whose trajectory 
lim^oo x(t) ends at x n . Note also that if point r G A, and the trajectory starting from r\ 
reaches r in a finite time t, then fj G A. This set defines a partitioning of space, where 
A n fl A m = when n ^ m and U ra y4 ra = f2. Finally, each basin A n is such that wherever 
the normal n to the bounding surface dA n is well-defined, n ■ Vp = 0. If p is the charge 
density, then A n are the Bader volumes; but this definition is applicable to any sufficiently 
smooth function with a discrete set of local maxima. As the definition of the basins A n 
derives from trajectories, it is not possible in general to determine if two neighboring points 
r and f* belong to the same or different basins based only on local information. 



FIG. 2. Schematic illustration of the weight method. The volume of the cell of a grid point flows 
to its neighbors with larger charge density magnitude. Flowing flux is shown as directional map, 
either flowing from X to X' as red arrows or flowing from X' to X as blue arrows. 

Figure [2] shows the reformulation for an approximate fractional partitioning of real- valued 
function evaluated at a set of discrete points, X. The grid points X partition space into 
Voronoi polyhedra 10 Vx covering each grid point X, where a point in space f belongs to the 
volume Vx if X is the closest point in Cartesian space to r. Each polyhedra is defined by the 
nearest neighboring points X' that are a distance lx-+x' away; the Voronoi polyhedron at X 
has facets dVx->x' with normal hx-^x' pointing from X to X' and area ax^x' ■ Moreover, 
the facet is at the midpoint between X and X' . Our goal is to define for each grid point 
X, a "weight" w A (X) between and 1 such that J2a wA (X) = 1 for all X, and the discrete 
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approximation to the integral over the basin A 



J A 



I rf 3 r/(r)^^^(X)/(X) 



(1) 



converges quadratically in the grid spacing for smooth functions f{r). The weight, in this 
case, is the fraction of points in Vx whose trajectory ends in the basin A. Note that if 
the points X form a regular periodic grid, the Voronoi volumes, facet areas, and neighbor 
distances need only be computed for the Wigner-Seitz cell around a grid point. 

To transition from the continuum definition of spatial partitioning to our Voronoi par- 
titioned definition, we introduce the continuum probability density for our trajectories, 
P(r, t). From the trajectory equation, the probability flux at any point and time is j(f, i) = 
P(r,t)Vp(f). Then, the probability distribution evolves in time according to a continuity 
equation 



This equation represents the combined evolution of a distribution of points in space; we use 
it to determine how the points in Vx distribute to neighboring volumes Vx 1 - Define the 
volume probability 



dP(f, t) 
dt 



+ V • (Pif.^Vpif)) = 0. 



(2) 




(3) 



then the evolution from the initial condition 





is given by 





Px(t)J2 T x^x> 



X' 
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where R(u) = u6{u) is the ramp function, so that tx^x> > and is zero when px> < px', 
this is a consequence of our initial conditions where P(r, 0) is only nonzero in the interior of 
Vx- The first approximation ignores spatial variation of P(r,t) through the volume Vx (an 
error linear in the grid spacing), and the second approximation ignores spatial variation of 
Vp along a facet dVx->x'i an d approximates the gradient at the midpoint between X and 
X' with the finite difference value (also with an error that is linear in the grid spacing). The 
solution to Eqn. [5] is Pxif) = Gxp(—tJ2x ,T x^x')- For that solution, the time-integrated 
flux of probability from Vx to Vx> through the facet dVx^x' is 



where we have used the same approximations as above. This flux defines the total fraction of 
points inside Vx that transition to volume Vx> through dVx->x'- Note that ^2 X ' Jx-+x' = 1 5 
unless X is a local (discrete) maxima, where px > Px> for all neighbors X' . Finally, as the 
weight w (X) represents the volume fraction of points in volume Vx whose trajectory ends 
inside basin A, then 



Note that if for all X' where p(X') > p(X), J2a wA ( x ') = x > then as Y,x> J x-+x' = 1, 
Eqn. [7] guarantees that ^2a wA ( x ) = !• Appendix lAl shows that the error in the weight of 
linear order in the grid spacing produces a quadratic order error in the integration. 

Forward substitution solves Eqn. [7| after the grid points are sorted from highest to lowest 
density p(X). Sequentially, each point X is either 

1. A local maxima: p(X) > p(X') for all neighbors X'. This grid point corresponds to a 
new basin A, and we assign w A (X) = 1. 

2. An interior point: for all X' where p{X') > p(X), the weights have been assigned 
and w A (X') = 1 for the same basin A. Then Eqn. [7] assigns X to basin A as well: 





J2x> a x^x't x \x>R(.Px' ~ Px) 




(7) 



X' 



w A (X) = 1. 
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3. A boundary point; with weights between and 1 for multiple basins assigned by Eqn. [7J 

Then w A (X) is known from w (X') where p(X') > p(X) for each basin A (as Jx^x' 7^ 
only if p(X') > p(X)). Note also that the weight for a particular basin A n is assigned without 
reference to any other basin A m \ once the set of time-integrated fluxes Jx-^x' are known 
and the densities sorted in descending order, the solution for each basin is straightforward, 
and Eqn. [7J is only needed on the boundary points. 

This algorithm solves several issues with the near-grid method. It requires no self- 
consistency, which improves the computational scaling. Moreover, the introduction of 
smooth functions that define the volume fraction of points in each basin produces less error 
and faster convergence with additional grid points. The algorithm is also readily applicable 
to non-uniform grids, such as an adaptive meshing scheme — it only requires computation of 
the Voronoi volumes and facets for the grid points. In one dimension, Eqn. [1] has quadratic 
convergence in the grid spacing (c.f. Appendix [AJ; we now demonstrate the quadratic 
convergence and improved integration accuracy for three dimensional problems. 

III. EVALUATION OF NUMERICAL CONVERGENCE 

One determination of the accuracy of Bader volume integration is the vanishing of the 
volume integration of the Laplacian of charge density V 2 p(r). The non-zero value of the 
Laplacian of charge density integration within each Bader volume is our atomic integration 
error, and can be used as an estimate of the error of the integration of the kinetic energy. We 
construct the zero flux surface of the gradients of charge density and evaluate the integration 
error with both the weight- and near-grid methods for several cases. First, we consider an 
analytic charge density with known boundaries in an orthogonal and a non-orthogonal cell. 
Next, we calculate real systems: an ionic compound, and a semiconductor. We also evaluate 
the Bader charge of Na atom in NaCl crystal by integrating the charge density within Bader 
volume, and compare the convergence with the near-grid method. 

A. Gaussian densities 

Figure [3] shows an example of misassignment of the grid points to basins from the near- 
grid method. Misassignment occurs for the grid points close to the dividing surfaces with 
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the gradients of charge density almost parallel to the surfaces. In this example, a three di- 
mensional model charge density is constructed from three Gaussian functions in simple cubic 
unit cell, p(r) = £\ =13 e^ r ~^ 2 / w \ The f { are (0.25iV, 0.25N, OAN), (0.5N,0.5N,0.5N), 
and (0.75N, 0.75iV, OAN), with width W = N/10. Figure [3] shows the charge density distri- 
bution on (110) plane. Due to the symmetry of charge density distribution, the true dividing 
surfaces along charge density saddle points are known analytically and shown as two black 
lines on (110) plane. The grid points marked by orange circles are assigned to the wrong 
basins by the near-grid method, different from the partition of spatial points by true divid- 
ing surfaces. The gradients of charge density shown in arrows for these misassigned surface 
points have small normal components, and we believe this is the cause of the misassignment. 

To test integration accuracy beyond simple cubic lattices, we map this model charge 
density onto a FCC unit cell shown in Figure HI The three dimensional model charge 
density is constructed from three Gaussian functions, p(r) = ^2 i=1 3 e(~ r ~ r ^ 2 / w2 . The are 
located at (0.25iV, 0.25iV, 0.4iV); (0.5N, 0.5N, 0.5N); (0.75N,0.75N,0AN) where N 3 is the 
number of grid points in the FCC unit cell and W = N/10. We vary N from 20 to 100. 
The Voronoi cell of FCC lattice has 12 neighbors, where all facets have the same area. The 
atomic weights on every grid represents the fraction of Voronoi volume of that grid point 
flowing to specific atom through its neighbors. By calculating on a set of grid sizes, one 
obtains the maximal atomic integration errors from the near-grid method and the weight 
method. 

Figure H] shows a reduction in error of three orders of magnitude from the near-grid 
method. Fitting data to a non-linear function y = aN~ r gives a convergence rate of 0.71 
for the weight method, and 0.45 for the near-grid method. The exponent of 0.71 is close to 
the 2/3 expected for quadratic convergence, and 0.45 is close to the 1/3 expected for linear 
convergence. The weight method has both better absolute error and converges faster than 
the near-grid method; in addition, there is no crossover point at large grid spacing where 
near-grid has smaller errors. 

B. Titania bulk 

For a real charge density, we perform DFT calculations on Ti02 bulk by use of the pro- 
jector augmented wave (PAW) 11 method, the GGA with PBE functional for the exchange- 
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FIG. 3. Charge density distribution constructed from three Gaussian functions, and basin identi- 
fication with near-grid and the weight method. The grid size is N = 60, and the charge density 
is shown on (110) plane; this is coplanar with the centers of the three Gaussian functions. The 
true dividing surfaces are indicated by two black lines due to symmetry. The basin assignment 
of grid points is given: red dots for ion I, green dots for ion II, and blue dots for ion III; for the 
weight method, a single color is assigned to the maximum weight at each point. Basin assignment 
from the near-grid method is given in the middle panel; basins with maximal weight on every grid 
points from the weight method are indicated in the right panel. Orange circles in the middle panel 
indicate the grid points misassigned by the near-grid method and corrected by the weight method. 
Arrows in the bottom panel denote the directions of the gradients of charge density, which can be 
used to verify the correctness of basin assignment. 

correlation energy. Density-functional theory calculations are performed with vasp^*^ us- 
ing a plane-wave basis with the projector augmented-wave (PAW) method,— with potentials 
generated by Kresse.— Atomic configurations for Ti and O are [Ne]3s 2 3p 6 4s 2 3<i 2 with cutoff 
radius 1.22A, and [He]2s 2 2p 4 with cutoff radius 0.58A, respectively. We use a plane-wave 

10 




FIG. 4. Comparison of the near-grid method and the weight method on atomic integration errors 
in a FCC cell. The maximal atomic volume integrations of the Laplacian of charge density within 
Bader volumes using the near-grid method and the weight method are denoted by squares and 
circles, respectively. We calculate charge density grids ranging from 20 3 points to 100 3 . Our algo- 
rithm gives atomic integration errors three orders of magnitude lower than the near-grid method, 
and converges faster than the near- grid method. 

basis set with cut-off energy of 900eV. The tetragonal unit cell of rutile Ti02 (see Fig- 
ure [5]) contains two Ti atoms and four O atoms. Monkhorst-Pack k-point method with 
4x4x6 k-points for six-atom cell is used for Brillouin-zone integration with a Gaus- 
sian smearing of O.leV for electronic occupancies. Theoretically optimized lattice constant 
are a = 4.649A, c = 2.970A, u = 0.305 agreeing with experimental lattice constants of 
a = 4.584A, c = 2.953A, u = 0.305^. A set of charge density grids ranging from 45 x 45 x 30, 
60 x 60 x 40, 75 x 75 x 50, 90 x 90 x 60, 120 x 120 x 80 points to 150 x 150 x 100 are calculated. 
For the energy cutoff of 900eV, a grid of 45 x 45 x 30 is required to eliminate wrap-around 
errors, and is the minimum size used by an accurate vasp calculation. 

Figure |5] shows maximal atomic integration errors as a function of grid sizes. The weight 
method gives maximal atomic integration error one order of magnitude lower than the 
near-grid method systematically. The atomic integration error larger than l.OeV on the 
minimal grid size 45 x 45 x 30 from the near-grid method is unacceptably large. Again, 
the convergence rate of the error goes as ~ 2/3 for the weight method — corresponding 
to quadratic convergence — and ~ 1/3 for the near-grid method — corresponding to linear 
convergence. Both the improved error and faster convergence allows for more accurate 
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density integration with fewer grid points than near-grid. 
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FIG. 5. Maximal atomic integration error on rutile T1O2 with respect to the charge density 
grids. A set of charger density grids ranging from 45 x 45 x 30 points to 150 x 150 x 100 points 
are calculated. The weight method reports maximal atomic integration error at least one order 
of magnitude smaller than the near-grid method. The weight method is practically useful for 
calculation of small grid size. 

C. NaCl crystal 

In this example, we evaluate the Bader charge (valence electron density integration within 
Bader volume) of Na atom in NaCl crystal by integrating the charge density within Bader 
volume, and compare the value with the near-grid method. We perform DFT calculations 
by use of the PAW method, the GGA with PW91 functional for the exchange-correlation 
energy. Atomic configurations for Na and CI are [He]2s 2 2p 6 3s 1 with cutoff radius 0.77 A, and 
[Ne]3s 2 3p 5 with cutoff radius 1.00 A, respectively. A plane-wave basis set with cut-off energy 
of 500eV is applied. The NaCl unit cell contains 4 Na atoms and 4 CI atoms. Monkhorst- 
Pack k-point method with 3x3x3 k-points for eight-atom cell is used for Brillouin-zone 
integration with a Gaussian smearing of 0.2eV for electronic occupancies. The optimized 
lattice constant of 5.67Aagrees with the experimental lattice constant of 5.64A. A set of 
charge density grids of 60 3 , 80 3 , 100 3 , 120 3 , to 180 3 points are calculated. 

Figure [6] shows the maximal atomic integration error as a function of various grid sizes. 
The weight method again shows maximal atomic integration error at least one order of 
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magnitude lower than the near-grid method systematically. The scaling of the error goes 
as the ~ 2/3 power for the weight method, showing continued quadratic convergence, while 
the near-grid method error scales as the ~ 1/3 power, which is linear convergence. 

Figure [7] shows that Bader charge of Na atom evaluated on various charge grids. The 
weight method computes a Bader charge of Na atom slightly larger than the near-grid 
method. Fitting the data to p = po + jM—, we find converged Bader charge values of 0.878e, 

grid 

0.881e, for the near-grid method and the weight method, respectively. We believe this is due 
to a systematic misassignment for the near-grid method, as shown for the Gaussian charge 
density case. This suggests that the misassignment may not be improved by increasing the 
density of grid points in the near-grid method. This suggests that a "divide and conquer" 
approach using continually refined grids can face potential difficulty For 60 3 grid points, 
the near-grid method underestimates the Bader charge by O.Ole, while the weight method 
underestimates it by 0.005e, again showing faster convergence. 




Grid Points 



FIG. 6. Comparison of the near-grid method and the weight method for maximal atomic integration 
error of NaCl crystal. A set of charge density grids ranging from 60 3 points to 180 3 points are 
calculated. Comparing to the near-grid method, the weight method reduce the integration error 
remarkably. 



IV. COMPUTATIONAL EFFORT 

The weight method is computationally efficient, requiring overall effort that scales linearly 
with the number of grid points. The total computer time is comprised of two primary 
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FIG. 7. Comparison of the near- grid method and the weight method on convergence of Bader 
charges of Na in NaCl crystal. The Bader charge of Na is calculated for a set of density grids 
ranging from 60 3 points to 180 3 points. Both methods give monotonic, and smooth convergence. 



tasks: the sorting of charge density costs O(NlogN) with N grid points, and the atomic 
weight evaluation on the sorted grid points beginning from grid point with maximum density 
requires at most iV x iV atom computer time. The computational effort is smaller than that, as 
only the surface grid points which have fractional atomic weights require iV atom calculations, 
while each interior grid point require only one calculation. Generally, the number of surface 
grid points is a small fraction of the number of total grid points, and scales as N 2 ^ 3 . For 
example, the ratio of the number of surface grid points to the number of total grid points is 
14% in NaCl crystal with total grid sizes 60 3 . 

In a calculation with charge density grid sizes approaching 10 7 -10 8 grid points and up 
to hundreds of atoms in large supercells, we find our algorithm is not only more accurate, 
but more efficient than the near-grid method. Both methods scale linearly with the number 
of grid points. Figure M shows the linear scaling of computer time required to analyze the 
charge density grid for an eight-atom NaCl with the number of grid points. The improved 
efficiency of our algorithm appears to originate from the lack of a self-consistent refinement of 
basin assignment. Comparing to the near-grid method, which needs refinement integration, 
our weight method has small prefactor, although both are linearly scaled. 
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FIG. 8. Computer time required to analyze the charge density grid for an eight-atom NaCl cell. 
The calculations were performed using an Intel Core2 Quad CPU Q6600, with a clockspeed of 
2.40GHz. The computer time scales linearly with respect to the number of charge density grid 
sizes with the weight method, as with the near-grid method. The weight method has a smaller 
prefactor than the near- grid method. 

V. CONCLUSIONS 

We develop a weight method to integrate functions defined on a discrete grid over basins of 
attraction (such as Bader volumes) in an efficient and accurate manner. The weight method 
works with the density on a discrete grid and assigns volume fractions of the Voronoi cell 
of each grid point to surrounding basins. Starting from the local density maxima, all the 
grid points are sorted in density descending order. Grid points can then be fractionally 
weighted from the weights of its neighbors with larger density. This method depends upon 
the formulation of flow across that dividing surfaces between the cells of two neighboring 
grid points, and can be applied to uniform or non-uniform grids. 

We perform tests on model three-dimensional charge density constructed from Gaussian 
functions in FCC cell. The weight method shows that the atomic integration error is in- 
versely proportional to the 2/3 power of grid points, while the integration error is inversely 
proportional to the 1/3 power of grid points using the near-grid method. We also perform 
tests on more realistic systems, such as Ti02 bulk and NaCl crystal. In both cases, the 
weight method reports maximal atomic integration error at least one order of magnitude 
smaller than the near-grid method systematically. Furthermore, we calculate the Bader 
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charge of NaCl crystal using these two methods, both give monotonic, and smooth conver- 
gence with respect to the increasing grid sizes, while they converge to slight different values, 
by 0.003 e. The weight method is more accurate than the near-grid method that require 
very fine grids. 
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Appendix A: Quadratic error in one-dimension 

The weight method for integration of the Bader charge volume has error that is quadratic 
in the grid spacing in one dimension. Consider the charge density p(x) evaluated on a regular 
grid with spacing h. In Eqn. [TJ there are only two grid points where w A (X) is not exactly 
or 1; these are the boundary points, and each is adjacent to a fixed point. In one dimension, 
the contributions to Eqn. [TJthat could produce errors linear in h only come from those points; 
the integration of the interior produces a total error that is quadratic in h. Hence, without 
loss of generality, we consider a single boundary point, and show that its contribution to 
Eqn. [1] produces an error that is of the order h 2 , rather than h. 

Let X be an boundary point, where the basin A lies to its left. This requires that 

p(X - h) > p(X), and p(X + 2h) > p(X + h). Finally, in order for w A (X) to not be 

identically 1, p(X) < p(X + h). This means that there is a point X + 5 for 5 G [0, h] such 

that p'(X + 5) = 0. Then, the flux from Eqn. [6] is 

p(X -h)- p(X) 
X ^ X - h ~ p(X -h) + p(X + h)- 2p(X) [Ai) 

and Jx^x+h = l-Jx^x-h] finally, as p(X) < p(X+h), J x +h^x = 0. Then, w A (X-h) = 1, 
w A (X + h) = 0, and so w A (X) = Jx^x-h- Finally, the contribution to Eqn. CD from X is 

rX+5 

/ f{x)dx « hJ x ^ x _ h f(X) (A2) 

JX-h/2 
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To evaluate the integration error, we use a Taylor expansion for p and / around the grid 
point X. We write p (n) = d n p/dx n (X) and f {n) = d n f/dx n (X). Note that pW has to scale 
as h in order for the dividing point X + 5 to lie between X and X + h. The Taylor expansion 
of w A (X) = Jx^x-h from Eqn. lAll to linear order in h is 

"1 pWh- r 



w A (X) ft 

so our integration contribution is 

"1 pW/i- 1 



(2) 



o(3) 

h-^ + 0{h 2 ) 



6p(2) 



(2) 



f{X)-h 2 ^f(X) + 0(h» 



(3) 



(A3) 



(A4) 



To find the true value of the expression, we need to determine 5 to at least quadratic order 
in h: write 5 = S^h + S^h 2 , and we have 



p '(X + 5) =p« + 5 ■ p< 2 ) + ^ 2 • p( 3 ) + 0(h 3 ) 
=h<jP>hr x ) + M«p( 2 ) 

+ ^ 5 (2) p (2) + 1^2 ( 5 (D) 2 p (3) +0(/r 3 

2 



(A5) 



which is solved by 



5 {i) 



2(p( 2 )) 3 

With our quadratic approximation for 5, we can integrate f(x) as 



(A6) 



X+S 



X-h/2 



f(x)dx =(x- X)f(X) + -(x- Xff^ + 0({x - Xf) 



x+s 



X-h/2 



-h 



P 



P 



(2) 



+ 2 V P (2) 



/(*) 

2 

7d) 



(A7) 



p (3) 



+ 0(/i 3 



which agrees with the contribution from our weight integration in Eqn. 
of order h 2 . As a special case, consider f(x) = 1/h; then the integral 



up to an error 



x+s 



dx = w A {X) + 0(h) (A8) 

I X-h/2 

which shows that the weight is the volume fraction of the Voronoi volume belonging to basin 
A to first order in h. 
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